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Abstract 

In this survey, we give a friendly introduction from a graph theory perspective to the g-state 
Potts model, an important statistical mechanics tool for analyzing complex systems in which 
nearest neighbor interactions determine the aggregate behavior of the system. We present the 
surprising equivalence of the Potts model partition function and one of the most renowned 
graph invariants, the Tutte polynomial, a relationship that has resulted in a remarkable synergy 
between the two fields of study. We highlight some of these interconnections, such as com- 
putational complexity results that have alternated between the two fields. The Potts model 
captures the effect of temperature on the system and plays an important role in the study of 
thermodynamic phase transitions. We discuss the equivalence of the chromatic polynomial and 
the zero-temperature antiferromagnetic partition function, and how this has led to the study 
of the complex zeros of these functions. We also briefly describe Monte Carlo simulations com- 
monly used for Potts model analysis of complex systems. The Potts model has applications as 
widely varied as magnetism, tumor migration, foam behaviors, and social demographics, and 
we provide a sampling of these that also demonstrates some variations of the Potts model. We 
conclude with some current areas of investigation that emphasize graph theoretic approaches. 

Keywords: Statistical Mechanics, Tuttc Polynomial, Potts Model, Ising Model, Monte 
Carlo Simulation, Chromatic Polynomial. 



1 Introduction 



The Potts model of statistical mechanics models how micro-scale nearest neighbor energy interac- 
tions in a complex system determine the macro-scale behavior of the system. This model plays an 
important role in the theory of phase transitions and critical phenomena in physics, and has ap- 
plications as widely varied as adsorption of gases on substrates, tumor migration, foam behaviors, 
and social demographics. If we generalize the regular lattice (on which physicists normally consider 



the Potts model) to an abstract graph, the g-state Potts model partition function is an evaluation 
of one of the most renowned graph invariants, the Tutte polynomial [Tut47, 53, 67, 84]. The Potts 
model is in fact equivalent to the Tutte polynomial if both q and the temperature are viewed as 
indeterminate variables. 

Here, we give a friendly introduction to the interconnections between the g-state Potts model 
partition function of statistical mechanics and the Tutte and chromatic polynomials of graph theory. 
In some respects, this paper complements the excellent survey of Welsh and Merino [WMOO]. 
Where [WMOO] is directed toward the physicist familiar with the Potts model who desires an 
introduction to the Tutte polynomial and its properties, here we hope to engage the graph theorist 
with an accessible introduction to the Potts model. Ideally, this paper will generate further interest, 
particularly from the graph theoretical perspective that has proven so productive in this rapidly 
developing and important area. 

The Potts model, from the mid- 1900s, builds on the seminal work of Ernst Ising [Isi25]. The 
Ising model of magnetic behavior features nearest neighbor interactions between spins at each point 
on a lattice, where the spins can assume either of two values corresponding to magnetic polarization. 
From these local interactions, the aggregate global properties of the system can be studied. Central 
among these are phase transitions, that is, critical temperatures around which a small change in 
temperature results in an abrupt change in the magnetism of the system. For the Ising model, 
there is no phase transition in the one dimensional case, and the transition has been determined 
exactly for the two dimensional square lattice. 

Intrigued by a related model due to Ashkin and Teller [AT43] , Cyril Domb suggested the study of 
what is now called the q^-state Potts model to his Ph.D. student, Renfrey B. Potts, who developed 
the beginnings of the theory in his 1952 doctoral thesis [Pot52]. The Potts model generalizes 
the Ising model by allowing q different spin values. Important thermodynamic functions such 
as internal energy, specific heat, entropy, and free energy, may be derived from the Potts model 
partition function. From a physics standpoint, one of the main reasons for the strong interest in 
the Potts model is that, for q = 3 and g = 4, it exhibits a continuous phase transition between 
high- and low-temperature phases with critical singularities in thermodynamic functions different 
from those of the Ising model. However, since its inception, myriad applications of the Potts model 
have emerged, and its usage now spans all areas of the sciences. 

We present the essential concepts in this area, addressing the natural questions of: What is 
the Potts model and its partition function? How are the Potts model partition function and the 
Tutte polynomial related? How does the Potts model capture phase transitions in thermodynamic 
functions? What is the relationship between the Potts model and the chromatic polynomial? What 
is the computational complexity of the partition function? How are Monte Carlo simulations used 
for the Potts model? Why is this model generating so much current interest? And finally, what 
are some current research directions in this field that emphasize graph theoretical approaches? 



2 The g-state Potts model 

Let G be a graph and 5* be a set of q elements, called spins. In the abstract, the spins may be 
numbers or colors, but typically they are values relevant to some specific application. For example, 
in studying uniaxial magnetic materials, q = 2, and the possible spins are +1 and —1. In a foam 



model, there may be thousands of spins, one for each bubble in the foam. In many applications, 
the graph G is taken to be a regular lattice, but this assumption is not necessary. 

A state of a graph G is an assignment of a single spin to each vertex of the graph. The 
Hamiltonian is a measure of the energy of a state. We begin with the simplest formulations, where 
the interaction energy (which may be thought of as a weight on each edge of the graph) is a 
constant J, and the Hamiltonian depends only on the nearest neighbor interactions (without any 
external field or other modifying forces). The model is called ferromagnetic if J is positive and 
antiferromagnetic if J is negative. If J is positive (respectively negative), the spin-spin interaction 
favors equal (respectively unequal) values of the spins on adjacent vertices. 

We will see shortly that the two Hamiltonians in the following definition generate essentially 
equivalent formulations of the Potts model partition function. 

Definition 2.1. Two common formulations of the Hamiltonian: 



ij&E[G) 



and /i2(a;) = J ^ (1 - 5 ((Tj, cjj)), 
ij(iE(G) 



where a; is a state of a graph G, Uj is the spin at vertex i, and 5 is the Kronecker delta function. 

For example. Figure [1] gives a graph state uj of the 4x4 square lattice with two choices of spin 
(black or white) at each vertex, with hi{uj) = — IIJ and h2{u>) = 13 J. Note that, up to the minus 
sign, hi counts the edges with the same spins on their endpoints, and /12 counts the edges with 
different spins on their endpoints, so, for any state w, /i2 (to) = J \E {G)\ + hi (uj). 
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Figure 1: A State of the 4x4 Square Lattice 

Definition 2.2. The Potts model partition function: Given a set of q spins and a Hamiltonian 
hi for i = 1 or 2, the g-state Potts model partition function is Pi (G) = ^ exp(— /3(/ii(a;))). Here 
the sum is over all possible states to of G, and P = k/T, where T represents the temperature of the 
system, and k = 1.38 x 10~^^ joules/Kelvin is the Boltzmann constant. 

The temperature T = 1/(k/3) is an important variable in the model, although it need not 
represent physical temperature, but may be some other measure of agitation or volatility relevant 
to the particular application (economic factors in a sociological model for example). Also, the 
product J (3 occurs often and customarily appears in the physics literature as Jf3 = K; we will 
adopt this convention when convenient. 

The Potts model partition function is the normalization factor for the Boltzmann probability 
distribution. For systems such as the Potts model that follow Bolzmann distribution laws, the 
number of states with a given energy (Hamiltonian value) are exponentially distributed. Thus, 



the probability of the system being in a particular state w at temperature T is: Pr (zu, /3) = 
exp(-/3/ij(ro))/ ^ exp(-/?/ii(a;)). 

Since the two different formulations of the Potts model partition function are each natural to 
use in different contexts, the following observation that one is simply a scalar multiple of the other 
facilitates translating theoretical results from one context to the other. 

Observation 2.3. 

^2(G';g,/3)=X:exp(-/3/i2M) = 

E exp (-/? ( J |£; (G)| + /ii H)) = exp (-i^ |ii; (G)|) Pi (G; (/, /3) . 

As a quick example, consider a single square, with two possible spins (white and black) at each 
vertex. The possible states (up to rotation) and their Hamiltonians using h\ are shown in Figure [2j 
Thus, the partition function for this graph is P\ (G) = 12exp {2K) + 2exp {4:K) + 2. 
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Figure 2: Possible States for C4 

As noted following Definition 2.2, the probability of a particular state w occurring is exp(— /9/ij(ci7))/Pj(G). 
Most importantly, since [3 = k/T, this probability is temperature dependent. Continuing the 
example in Figure [21 we take K = l/T, so J = k is positive, and consider the probability 
of the all-black state occurring at different temperatures. Since J is positive, and we are us- 
ing /ii, the all-black state is one of the two lowest energy states. Thus, we would expect the 
system to favor the two low energy states equally at low temperatures, and be equally likely 
to be in any of the sixteen possible states at high temperatures. If we let ojh be the all-black 
state, then the probability of the system being in the all-black state, as a function of tem- 
perature, is Pr(u;6,T) = exp(4i^)/(12exp(2i(') -|- 2exp(4i^) -|- 2). Evaluating this at various 
temperatures illustrates how the model captures the expected temperature dependent behavior: 
Pr(a;6, IQ-^) = 1/2 = 0.50, Pr(a;fe, 2.29) = 0.19; Pr(u;b, 10^) = 1/16 = 0.0625. 



3 Relating the Potts model and the Tutte polynomial 

From a graph theory perspective, one of the most remarkable aspects of the Potts model is its 
connection with the Tutte polynomial, one of the best known graph invariants. See Fortuin and 
Kasteleyn [FK72] for the nascent stages of this discovery, and Wu [Wu82, Wu88] for further discus- 
sion from a physics viewpoint. More recent mathematical physics studies, such as Shrock [ShrOO], 
Sokal [SokOO], and Welsh and Merino [WMOO], give further exposition, while some relevant math- 
ematical reviews include Tutte [Tut84], Biggs [Big93], Bollobas [Bol98], and Welsh [Wel93]. The 
Tutte polynomial [Tut47, 53, 67] has a rich history and a wide range of applications. We mention 
here only a very few properties necessary to establish the relation between it and the Potts model. 



and refer the reader to the relevant chapters of Welsh [Wel93] and Bollobas [Bol98] and also to Bry- 
lawki [BrySO] and Brylawski and Oxley [B092] for an in-depth treatment of the Tutte polynomial, 
including generalizations to matroids. 

The Tutte polynomial, t{G;x,y), is a two- variable graph invariant that may be defined by a 
linear recursion relation in terms of deleting and contracting edges. Recall that an edge is deleted 
from a graph G by removing the edge, but not its incident vertices. This is denoted by G — e. An 
edge is contracted in G by removing the edge and coalescing its incident vertices. This is denoted 
by G/e. An isthmus (a.k.a. a bridge, cut-edge, or co-loop) is an edge whose deletion increases the 
number of components of G. A loop is an edge where both ends of the edge are incident with the 
same vertex. 

Definition 3.1. Deletion-contraction definition of the Tutte polynomial: t{G;x,y) = t{G — 
e; X, y) + t{G/e; x, y) if e is neither an isthmus nor a loop, and t {G\ x, y) = x^y^ if G consists of i 
isthmuses and j loops. 

As a quick example, in Figure [3]we calculate that t (C4; x,y) = x'^ -\-x'^ + x + y, where C4 is the 
cycle on four vertices (a square). In the diagram, the labeled edge is deleted and contracted in the 
next step, and a graph consisting of only isthmuses and loops is evaluated as a monomial in x and 
y in the following step. 



./.Q 



X + X ■*- U ^ ^-# = X ^ + X + X + y 



Figure 3: Calculation for C4 

A surprising (and essential) property of the Tutte polynomial is that the result of the recursion 
process is independent of the order in which the edges are contracted and deleted. Recall that for 
a graph G = {V,E), a spanning subgraph is G' = {V,A) where A is a subset of E. Also, k{G) 
is the number of connected components in G and c(G) is the number of independent cycles in G. 
One way to prove that t{G; x, y) is independent of the order of contraction and deletion (and hence 
well-defined as given in Definition 3.1) is to express it as a sum of terms, each of which is evaluated 
for a given spanning subgraph. This is analogous to the cluster representation for the Potts model 
in Fortuin and Kasteleyn [FK72]. In particular, t{G; x,y) = Y^{x — 1)'=('^ )-KG)[y — ly^^ ), where 

G' 
the sum is over the spanning subgraphs G' of G (see Bollobas [Bol98], and Biggs [Big93]). Since 

the set of spanning subgraphs is unique, this proves that t(G; x, y) is well defined for a given 

graph G. Alternatively, using induction on the number of edges, one can show that t{G] x, y) is 

equivalent, up to a prefactor, to the dichromatic polynomial, defined by Tutte [Tut67, 47,53] as 

Z{G-u,v)= Z u'=(^^i^l^l~l^(^)l+'^(^). Specifically, 
ACE(G) 

^kiG)^\v\-kiG)^ fG;^,v + l]=Z {G; u, v) . (1) 



Theorem 3.1 gives another critical property of the Tutte polynomial. The Tutte polynomial is 
universal in that essentially any other multiplicative graph invariant that has a deletion-contraction 
reduction must be an evaluation of it. 

Theorem 3.1. If f{G) is a function on graphs with a, b having ah 7^ 0, such that 

A. f{G) = 1 if G consists of only one vertex and no edges, 

B. f{G) = af{G — e) + bf{G/e) whenever e is neither a loop nor a bridge, 

C. f{GH) = f{G)f{H) where either GH is the disjoint union of G and H or G and H share at 

most one vertex, 

then f is an evaluation of the Tutte polynomial of the form f{G) = a"'b'^t{G; xq /b, yo/ a) . Here 
n is the nullity of G, that is, n = \E{G)\ — \V{G)\ + k{G), and r is the rank of G, that is, 
r = \V{G)\ — k{G). Also, xq = /(i^2) and yo = f{L), where K2 is the complete graph on two 
vertices, and L is a single vertex with a single loop edge. 

This universality provides the basis for a proof of the connection between the Tutte polynomial 
and the Potts model partition function. If we consider the hi Hamiltonian in the Potts model 
intuitively, we note that, in a given state of the graph, if the endpoints of an edge have different 
spins, then the Kronecker delta value is zero, and the edge contributes nothing, so it might as 
well be deleted. On the other hand, if the endpoints have the same spin, they interact with the 
neighboring points in exactly the same way, so they may be coalesced, with the edge contracted. 
However, this edge does contribute to the Hamiltonian, so there is a weighting factor when the 
edge is contracted. 

Although the Potts model partition function seems likely to satisfy the conditions of Theorem 
3.1, and thus be an evaluation of the Tutte polynomial, it does not satisfy condition A: if there 
are q spins, the Potts model partition function of a single vertex is q, not 1. However, there is a 
very common device for applying Theorem 3.1 in such situations, namely introducing a factor of 
some term raised to the power of k (G), a trick employed in the proof of the following theorem. 
Although the connection between the Potts model partition function and the Tutte polynomial was 
first recognized by Fortuin and Kastelyn [FK72], the proof below is modeled on those in Welsh 
[Wel93] and Bollobas [Bol98]. 

Theorem 3.2 Let P {G;q, /3) = q-^^^^Pi{G;q,(3). 

Then P (G; q, (3) = 7;|V{G)|-fc{G)^ ^q. (^ ^ ^y^^ ^ ^ y), 

and thus Pi {G;q,(3) = qk(G)^\V(G)\-k{G)^ ^q. (^ ^ v)/v,v + 1), where v = e^^{J[3) - 1. 

Proof: For ease of reading, we will suppress the subscript 1 on P and h in the following proof. 
The proof consists simply of verifying that P satisfies the conditions of Theorem 3.1. P clearly 
satisfies condition A. 

For condition B, let e = {c, d} be an edge of G which is neither a loop nor an isthmus, and 
write s (c) and s (d) for the spins at c and d respectively. Then 



PiG;q,P) = q-'^^^ ^ exp {- f3h (uj)) 



a)£ states of G 
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ui£ states of G 
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exp (— /?/l (lo)) 
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-k(G) 



J 
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Y, exp(-/3/iH) 



ljS states of G 
y with s{c)=s{d) 



(2) 



/ 



Note that if log is any state of G, then there is a unique state wc-e of G — e where each vertex 
has the same spin as it has in ug- Also, if ug is any state of G — e with s (c) = s (d), then there 
is a unique state cog/c of G/e, where the vertex resulting from identifying c and d has the common 
value s (c) = s{d), and each other vertex has the same spin as it does in ujg-c- Furthermore, if 
s (c) 7^ s{d), then hiuJG) = h{uJG-e)', and if s (c) = s{d), then /i (wg) = h{ujG-e) + J- Thus, 
equation (2) becomes 



/ 



\ 



-KG) 



Y^ exp{-(3h{u;)) 



uj£ states of G—e 
\with s{c)j^s{d) 



I 



\ 



+ g-^(^) exp (J/3) 



Y exp(-/3/iH) 



/ 



Va;£ states of G—e 
with s(c)=s(d) 



/ 



The right-most term is nearly P {G — e;q,P), since e being neither a bridge nor a loop means 
that k (G) = k{G — e), but we are missing the states of G — e where s (c) = s (d). So we simply 
add and subtract them, getting 



/ 



-k{G-e) 



Y exp(-/3/iH) 



\ 



VuiG states of G—e 
with s(c)^s(d) 



+ q' 



-k(G-e) 



exp (J/3) 



/ 



/ 



VcjG states of G—e 
with s(c)=s{d) 



) 
I 



-k{G-e) 



( 



VtJG states of G—e 
with s{c)=s(d) 



\ 



Y exp(-/3/i(a;)) 

3f 
/ 



/ 



+ 



-fc(G-e) 



a;£ states of G—e 
\with s{c)=s{d) 



-k(G-e) 



(exp (J/3) - 1) 



\ 



Y exp(-/3/i(^)) 



tjg states of G—e 
\with s(c)=(d) 



+'?■ 



-fc{G-e) 



/ 



Y exp(-/3/i(a;)) 

\ 

y 

( Y exp(-/3/i(^)) 

\a;G states of G—e / 



Y exp{-(3h{Lo)) 



The second term is now P {G — e; q,(3). For the first term, note that since e is neither a bridge 
nor a loop, k{G — e) = k (G/e). Also, the states of G — e with s{c) = s (d) correspond exactly to 
the states of G/e, and furthermore a state of G — e with s (c) = s (d) has the same Hamiltonian as 
the corresponding state of G/e. Thus, the first term becomes 



-k{G/e) 



(exp(ir)-l) Y exp (-/3/i (a;)) = (exp (ir)-l)P (G/e; (?,/?) 



u>£ states of G/e 



Thus, if e is neither a bridge nor a loop, then P{G;q,(3) = P {G — e\q,(3) -\- (exp(i^) — 1) 
P (G/e; g, /?), which satisfies condition B with a = 1 and h = exp {K) — 1. 



For condition C we write G U H when G and H are disjoint, and G * H when G and H share 
a single vertex. Condition C is easily satisfied when G and H are disjoint, since in this case a 
state ujgvjh of G U H is just two independent states coq and ujh of G and H respectively, and 
h {ujgvjh) = h {ujg) + h {ujh)- With this, and noting that A; (G U H) = k{G) + k{H), it follows that 



P{GUH;q,P)=q 



~k{GUH) 



^ exp(-/3/ii(w)) 



states iij 
oiGUH 



-k(G)-k{H) 



J2 exp{-P{h{ujG) + h{ujH))) 

UJQ of 
Uj£j of 

-^(^) J2 exp (-/?/! (^g))<Z-'^^) E exp(-/3/i(a;H)) 



states uiQ of G and 
states uj[j of // 



states uJG of Gf states ljj:/ of H 

P{G;q,(3)P{H;q,(3). 

In the case where G and // share a single vertex n, a state ujg*h of G * H corresponds to two 
states ujg and ujjj which both have the same spin at u. Here, k{G * H) = k{G) + k{H) — 1. Note 
that the number of states of H in which s{u) = a is equal to the number of states in which s{u) = b, 
for any other spin b, since we can simply exchange the roles of a and b in any state. Thus, 



E exp{-l3h{uJH)) 

wh £ states of H 



q E ex.p{-f3h{ujH))- 



UJH S states of H 
with s{u)=a 



With this we have: 



P{G*H;q,l3)=q 



-k{G*H) 



Y, exp(-/3/iH) 



states u] 
olG*H 



q-HG)-k(H)+i J2 exp{-(3h{ujG))ew{-ph{uJH)) 



ojh G states of H 
with s{u) in ijjj{ 
equal to s{u) in loq 



-k{G)-k{H)+l 



E 



ojqG states of G 



\ 



X! exp(-/3/i(wG))exp(-/3/i((jH)) 



i^H S states of H 

with s{u) in cj^f 

y equal to s(u) in tjQ 



/ 



-k{G)-k{H)+l 



q ■--"' ■-■^"'•^ ^ g 1 ^ exp(-/?/i(a;G))exp(-/3/i((jH)) 

ti^G G states of G i^hS: states of H 

-k[G)-kiH) Y^ exp{-Ph{uG)) Yl exp{-Ph{uH)) = 

iL)G S states of G u)jj G states of H 

P{G;q,p)P{H;q,P). 



Thus, condition C is satisfied, and Theorem 3.1 apphes to P. To apply the conclusion of 
Theorem 3.1, it remains to compute P at a single isthmus and a single loop, that is, find P (K2) 
and P(L). 

For a loop, note that there are q states, and since both endpoints of a loop necessarily have the 
same value, the Hamiltonian of each state is 1. Thus, 

P{L;q,P)=q-' ^ exp (-/3 (- J • 1)) = g-^gexp (iT) = exp (A^ . 

q states 

For K2, there are q states where the spins on the endpoints are equal, giving a Hamiltonian of 1. 
Then there are q{q — 1) states where the spins on the endpoints are different, giving a Hamiltonian 
of 0. Thus, 

P [B- q, (3) = g-i {q {q - 1) exp {-K ■ 0) + q exp {K • 1)) = (exp (A') + q - 1) . 

We now apply Theorem 3.1 with a = 1, b = (exp(A') — 1), xq = (exp (A') + g — 1), and 
yo = exp (AT). To simplify the expression, we first set u = exp (A') — 1, so a = 1, b = v, yo = v + 1, 
and xq = q + V. Now Theorem 3.1 yields 






Thus, since P (G; q, (3) = q-^^^^Pi (G; q, (3), it follows that 

A (G;g,/3) = q>^(0)^\no)\-HG)^ j^^. i±^,. + l) . /// 

As an example, we recall from Figure[2]that Pi (G) = 12 exp(2A:)+2 exp(4A)+2, so Pi {K2, 2, (3) = 
12 (-y + 1) + 2 (f + 1) + 2, since there were q = 2 spins in the example. Also recall from Figure [3] 
that t{Ci]x,y) = x^ + x'^ + x + y. A quick calculation verifies that 2^v'^~^t (G; "^^^v + l) = 
12(w + l)V2(?; + l)V2. 

We then have the following two immediate corollaries, with Corollary 3.3 following from Obser- 
vation 2.3, and Corollary 3.4 from equation (1). 

Corollary 3.3. 

P2 (G; q, 13) = g'^W {v + l)-l^(^)l J^(^)l-'=(G)t (g; ^, exp (K) 



Corollary 3.4. 



Pi{G;q,P) = Z{G;q,v)= J^ q 



HA)JA\ 



v 



ACE(G) 



Thus, if q and v are viewed as indeterminates, the Potts model partition function is exactly 
equal to the dichromatic polynomial Z {G;q,v), and in fact the partition function is typically 
denoted by Z in the physics literature. Corollary 3.4 also leads to the property (not apparent from 

9 



the original definition) that the Potts model partition function is a two-variable polynomial with 
maximal degree in q equal to the number of vertices of G and maximal degree in v equal to the 
number of edges of G. Theorem 3.2 and Corollary 3.4 also show how the Tutte polynomial may 
be thought of as an analytic continuation of the Potts model with its positive integer values for q. 
The random cluster model of Fortuin and Kasteleyn [FK72] also extends the Potts model in this 
way to M"*", and is likewise an evaluation of the Tutte polynomial. 

One common extension of the Potts model, dating to the earliest work in the quantum theory 

of magnetism in the 1920's and 1930's, involves allowing the interaction energy to depend on the 

specific edge, rather than be constant throughout. With this, the Hamiltonian becomes h[Lo) = 

^ Jij5{(Ti,aj), where Jij (or Jg) is the interaction energy on the edge e = {i, j}. The partition 
{ij}&E{G) 

function then becomes P{G) = Z{G;q,v) = J2 q^^ ' Yl ^e, where Ve = exp(/3Je) — 1. See, for 

ACE{G) eeA 

example, Fortuin and Kasteleyn [FK72] and Baxter [Bax82] for an edge weight generalization of 
the Potts model, and also recent work by Sokal [SokOO, Ola]. 

In recent years, the Tutte polynomial has also been extended to incorporate edge weights. 
Here, these weights may also depend on whether the edge is contracted, deleted, or evaluated as 
an isthmus or a loop as the polynomial is recursively computed. In the most general case, however, 
care must be taken with a set of relations on small graphs (two edges in parallel, a cycle with 
three edges, and three edges in parallel) to assure that the resulting function is well defined. See 
Traldi [Tra89], Zaslavsky [Zas92], Bollobas and Riordan [BR99], and Ellis-Monaghan and Traldi 
[E-MT06]. The generalized partition function given above satisfies these relations, however, and 
thus the connection between the Potts model partition function and the Tutte polynomial extends 
to systems with edge-dependent interaction energies. 



4 Thermodynamic functions and phase transitions 

An important goal of statistical mechanics is to determine phase transitions, that is, critical tem- 
peratures around which a small change in temperature results in an abrupt, nonanalytic change in 
various physical properties. Roughly speaking, this phase transition temperature separates the two 
phases of the system. For temperatures above this critical temperature, the system (in the absence 
of an external biasing field) exhibits no long-range order, e.g., no spontaneous magnetization in the 
case of a magnetic system, while for temperatures below the critical temperature it does exhibit 
such order. (See Stanley [Sta71], and Plischke and Bergesen [PB06].) In the following discussion, we 
elide technical caveats concerning such details as the choice and growth of the lattices, interchang- 
ing limit signs, boundary conditions, and convergence, etc., in order to provide a broad picture 
of the general principles. However, treating these technicalities carefully can present significant 
challenges in determining phase transitions for various applications. Bearing this in mind, when 
we speak of taking the thermodynamic (or infinite volume) limit below, we mean specifying an 
appropriate infinite family of graphs, such as square lattices, and taking the limit of an expression 
as the size of the graphs goes to infinity. 

In the Potts model, important thermodynamic functions such as internal energy, specific heat, 
entropy, and free energy (denoted U, C, S, and F, respectively) may all be derived from the 
partition function, Z = Z {G;q,v). For example, the internal energy, which is the sum of the 
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potential and kinetic energy and is defined as U = -^J2^ i^) ^^P {~P^ ('^))) ™&y be expressed by 
~ l^g ■ Specific heat, or the energy required to raise a unit amount of material one temperature 
increment, is C = gj?. Entropy, a measure of the randomness and disorder in a system, is 5* = 
— K/3 — 2i-^ + Kin (Z). Finally, the total free energy is F = U — TS = —kT In (Z). 

It is convenient to work with the dimensionless, reduced free energy / = —(3F, so the reduced 
free energy per unit volume (or in this context, per vertex) is 

f {G;q,f3) = ^^^^ln{P^{G;q,(3)) . 

For a fixed graph G, this is clearly an analytic function in both q and T. Any failures of 
analyticity can only occur in the infinite volume limit, that is 

f{T;q,(3)= hm / In (A (G„; g,/?)) , 

n^oo \V [Gn)\ 

where P is an appropriate infinite family of graphs. 

If the resulting limit of the reduced free energy is an analytic function of T (generally q is fixed), 
then the model has no phase transitions since U, C, F, and 5 will be analytic as well. This happens 
for example for finite temperature /3 G [0, oo) in the Ising model on the one-dimensional lattice. 
Values of T where analyticity fails are critical temperatures corresponding to phase transitions. 
Since phase transitions are manifest as failures of analyticity in the thermodynamic limit of the 
reduced free energy, the goal is either to determine these points of nonanalyticity or to establish 
analyticity in some region, and here the behavior of the partition function is the key. Furthermore, 
if we want to know how the thermodynamic functions behave near a critical temperature, again 
understanding the partition function is essential. 

If Tc is a critical temperature, and we write t = (T — Tc)/Tc, then the goal is to express any of 
these thermodynamic functions in the form g (r), where g (r) is roughly equal to c\t\^ near r = 
for some p called the critical exponent. The critical exponents fall into a set of discrete universality 
classes, where universality means that the values of the critical exponents are independent of 
parameters such as the interaction energy J and the choice of lattice (although not its dimension). 

Phase transitions are broadly classified as first-order (discontinuous) or second-order (contin- 
uous, but nonanalytic). Second-order phase transitions are further classified by the critical ex- 
ponents, since if the phase transitions are continuous, the leading singular behavior as r ^ of 
the thermodynamic quantities may normally be written in the form c |r|^, where p is a positive or 
negative power (there are also cases where the singularity is non-algebraic). See Stanley [Sta71], 
Fisher [Pis 74], and Plischke and Bergesen [PB06]. 

One of the important features of the two-dimensional ferromagnetic g-state Potts model is that, 
for the thermodynamic limit of regular two-dimensional lattice graphs, it provides, within one 
model, a set of several different universality classes associated with second-order phase transitions 
depending on q, in particular for q = 3 and q = 4:, which generalize the Ising q = 2 case. For q > 5, 
the phase transition of this two-dimensional ferromagnetic Potts model is first-order. Onsager 
[Ons44] calculated an exact closed-form expression for the free energy of the Ising model on a 
square lattice (in the absence of an external magnetic field) in 1944 (reviewed by McCoy and Wu in 
[MW73]). For values of g > 3 , the free energy of the g-state Potts model for arbitrary temperatures 
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(either with or without an external magnetic field) has never been calculated exactly for (the infinite 
volume limit of) a lattice of dimension two or more. 

Just as studying the generalization to complex variables can lead to greater insight into prop- 
erties of functions of real variables, the generalization of the variable K = J/{kT) from real values 
(positive for ferromagnetic, negative for antiferromagnetic) to complex values has proven quite in- 
formative. Indeed, this is necessary in order to study zeros of the partition function. For fixed q, 
the accumulation set of the zeros of the partition function in the infinite volume limit form curves 
in the complex plane referred to as the phase diagram. For example, Fisher [Fis65] showed that for 
the Ising model on the square lattice, as the number of vertices goes to infinity in the infinite volume 
limit, the zeros of the partition function asymptotically merge to form two circles in the complex 
plane, \v\ = \/2 and |w + 2| = \/2. (Earlier, Lee and Yang had studied the corresponding extension 
to complex numbers of the magnetic field [YL52, LY52].) This leads to the study of regions in the 
plane of the complex-temperature variable v which are analytic continuations of the physical high- 
and low-temperature phases of the ferromagnetic and antiferromagnetic models. Determinations of 
the complex-temperature phase diagram for the (7-state Potts model for arbitrary q were given for 
finite- width, infinite-length strips in Shrock [ShrOO], Chang and Shrock [CSOO, 01b], Chang, Salas 
and Shrock [CSS02], Chang, Jacobsen, Salas and Shrock [CJSS04], Jacobsen, Richard and Salas 
[JRS06], and the shape of these phase diagrams in general is discussed by Biggs [Big02b]. 

The determination of these phase diagrams requires a particularly interesting combination of 
methods from mathematical physics, graph theory, complex analysis, and algebraic geometry (since 
the phase boundaries are algebraic curves). Given the exact partition function as a function of 
both arbitrary q and v, it is possible to determine regions of analyticity in the complex q plane as 
a function of the temperature variable v for both the ferromagnetic and antiferromagnetic cases. 
In addition to determining the phase diagram in v for fixed q and in q for fixed v, one can also 
determine it when q and v satisfy a given functional relation (see Chang and Shrock [CS06]). 
Complex-temperature zeros of Potts models for fixed values of q (beyond the known Ising case) 
have also been calculated on finite patches of two-dimensional lattices (e.g. Martin and Maillard 
[MM86], Martin [Mar91], Chen, Hu and Wu [CHW96], Matveev and Shrock [MS96], Kim and 
Creswick [KCOl]). 

Additionally, Salas and Sokal [SS97] have shown that the antiferromagnetic Potts model with q 
spins has no finite temperature phase transition on lattices where each vertex has degree less than 
q/2. 

Quite a lot is known in the ferromagnetic case as opposed to the antiferromagnetic case. For 
example, the value of the critical temperature for the ferromagnetic Potts model on the infinte 
volume limit of the square lattice has been determined to be kTc = J/ln(l + ^/q), and the critical 
behavior of the two-dimensional Potts ferromagnetic model is known (see Baxter [Bax82] and Wu 
[Wu82]). Simulations on the square lattice agree with this formula for the critical temperature, and 
there have been a number of studies of the Potts model on various 2- and 3-dimensional lattices, 
giving valuable approximations with sufficient accuracy for relevant applications. Wu [Wu82, 84] 
and Salas and Sokal [SS97], for example, provide a survey of results and approximations. Further 
insight into the critical exponents has been gained from the use of conformal field theory (Cardy 
[Car87], Di Francesco, Mathieu and Senechal [DiFMS96]). 

Reviews of the Ising model include McCoy and Wu [MW73]. Cipra [Cip87] includes a highly 
accessible treatment of the 1-dimensional Ising model and the existence and nature of a phase 
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transition for the 2-diniensional Ising model. This is discussed in a more general setting in Stanley 

[Sta71]. 



5 Extremal temperatures and the chromatic polynomial 

In addition to critical temperatures in (0,oo), and the straight-forward case of /? = (i.e. infinite 
temperature), the extremal case of zero temperature, in both the ferromagnetic and antiferromag- 
netic models, is studied. Of particular interest, with respect to our theme of the interconnections 
between graph theory and statistical mechanics, is the equivalence of the zero temperature anti- 
ferromagnetic Potts model partition function and the chromatic polynomial. We will first briefly 
mention the infinite temperature and the ferromagnetic cases, and then focus attention on the 
chromatic polynomial and its zeros. 

In statistical physics, the free energy, F = U — TS, is minimized in a system in thermal 
equilibrium. Note that at low temperatures, the minimization of the free energy becomes equivalent 
to the minimization of the internal energy U. For positive J this means that all spins take on the 
same value, while for negative J it means that adjacent spins must have different values. As T 
increases, the minimization of the free energy increasingly means the maximization of the entropy, 
S. As T approaches infinity, the minimization of the internal energy plays a negligible role relative 
to the maximization of entropy. 

At infinite temperature, i.e., when i^ = 0, we have Z {G; q, v) = g'^C-')!, so / = Inq. It is then 
possible to calculate the high-temperature Taylor series expansions for thermodynamic quantities 
as powers oi v = exp{K) — 1, so around v = 0. These expansions are carried out around v = 
because it is possible to systematically generate (e.g. via graphical techniques as in Nagle [Nag71] 
and Kim and Enting [KE79]), higher order terms in the Taylor series expansion about this point. 

For the ferromagnetic Potts model at T = 0, the system exists in a completely ordered state, in 
which all spins have the same value, so that as T ^ 0, and hence X — > oo , the partition function 
Pi {G;q,P) — > q exp (K \E (G)\) . For a regular lattice whose vertices have uniform degree 6, the 
reduced free energy per unit volume is then |y(G)p Inq + ^K, so as the number of vertices goes 
to infinity, this gives the free energy per vertex as simply F = —((5/2) J. The low-temperature 
series expansion can then be expressed as a series in exp(— i^). Reasonably accurate values for 
both the critical temperature and the critical exponents can then be extracted from the Taylor 
series expansion. 

This now brings us to the zero-temperature antiferromagnetic case and the chromatic polyno- 
mial. We first recall that a proper coloring of a graph G is an assignment of a color to each vertex 
of G so that any two adjacent vertices receive different colors. The chromatic polynomial, C{G; x), 
is a graph invariant that, when evaluated at a non-negative integer x, gives the number of ways 
to properly color the graph G using x colors. Consider for a moment an edge e of a graph G. 
The number of ways to color G — e (where there are no restrictions on the colors assigned to the 
endpoints of e) is equal to the number of ways to color G (where the endpoints must have different 
colors) plus the number of ways to color G/e (where the endpoints, now coalesced, must have the 
same color). This means that the chromatic polynomial may be computed recursively as follows: 

G{G -e;x) = G{G; x) + C{G/e; x) or, C(G; x) = C{G -e;x)- G{G/e; x), 

and C{G;x) = x", if G has n vertices and no edges. 
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Thus, the chromatic polynomial has a contraction/deletion reduction and hence by Theorem 
3.1 must be an evaluation of the Tutte polynomial. An argument almost identical to the proof of 
Theorem 3.2 shows that 

C{G; x) = (-l)l^(«)l-^(«)a;'=(«)t(G; 1 - x, 0). 

Some early reviews on the chromatic polynomial include those of Read [Rea68], Read and Tutte 
[RT88], and Biggs [Big93]. An extensive bibliography is available in Chia [Chi97], while Thomassen 
[ThoOl] and Dong, Koh, and Teo [DKT05] give recent comprehensive treatments. 

The connection between the Potts model and the chromatic polynomial occurs in the zero- 
temperature (limit as T — > oo) antiferromagnetic model using the /ii(w) = — J Yl ^{^i-,(^j) 

ijeE{G) 
formulation of the Hamiltonian. Since J is negative in the antiferromagnetic model, minimal 
energy states are those that generate a maximum number of zeros in the summation, i.e. those in 
which every edge has endpoints with different spins. If we think of the spins as colors, a minimum 
energy state is then just a proper coloring of the graph. 

We give two different ways to understand the translation between the zero temperature Potts 
model and the chromatic polynomial. One approach is to compare C{G;x) = (— l)l^('^)l~'"AG')2.fc(G) 
t{G; 1-x, 0) with the result of Theorem 3.2 that Pi (G; q, (3) = qKG)y\V(G)\-k(G)^ ^q. (^ _^ ^y^^ v + 1), 
where v = exp(J/3) — 1. We note that these are the same function precisely when v = —1, that is, 
when /? = oo, which is exactly the zero-temperature model. Another way to view this connection 
is by considering the summands of Pi{G;q,(3) = ^ exp (/3J^ 5 (o-j, ctj)). As T ^ 0, and hence 
P — > oo, a summand is except precisely when ^ 6{cTi, ctj) = 0, in which case it is 1. Thus Pi{G) 
simply counts the number of proper colorings of G with q colors. 

In the special case of the Potts antiferromagnetic model at T = 0, where the Potts model 
partition function reduces to the chromatic polynomial, there are power series expansions of the 
various thermodynamic functions. When T = 0, the model will be in one of its possible ground 
states. Ground state entropy is a measure of the residual disorder in the system, and it can be 
nonzero for sufficiently large q. In the infinite volume limit, the ground state entropy per vertex of 
the Potts antiferromagnetic model becomes 



n^oo |K (G„)| 

This is related to the ground state degeneracy per vertex, W, according to S = kIuVI^. Two 
exact results are Lieb's calculation in [Lie67] giving W = (4/3)^' ^ for q = 3 on the square lattice 
(see also Baxter, Kelland and Wu [BKW76]) and Baxter's calculation of W for general q for the 
triangular lattice (see Baxter [Bax87], Wu [Wu82] and also Blote and Nightingale [BN82] and 
Baxter [Bax86]). More generally, some calculations of Tutte polynomials for recursive families of 
graphs have been carried out in Shrock [ShrOO], Chang and Shrock [CSOO, Ola, 06], Chang, Salas 
and Shrock [CSS02], Change, Jacobsen, Salas and Shrock [CJSS04], and Jacobsen, Richard and 
Salas [JRS06]. 

A significant body of work has emerged in recent years devoted to clearing regions of the com- 
plex plane (in particular regions containing intervals of the real axis) of roots of the chromatic 
polynomial. Results showing that certain intervals of the real axis and certain complex regions are 
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free of zeros of chromatic polynomials include those given by Woodall [Woo92], Jackson [Jac93], 
Shrock and Tsai [ST97a,c], Thomassen [Tho97], Brown [Bro98], Sokal [SokOlb], Procacci, Scop- 
pola, and Gerasimov [PSG03], Choe, Oxley, Sokal, and Wagner [COSW04], Borgs [Bor06], and 
Fernandez and Procacci [FP] . One particular question concerns the maximum magnitude of a zero 
of a chromatic polynomial and of zeros comprising region boundaries in the q plane as the number 
of vertices \V\ — > cxd. An upper bound is given in [SokOlb], depending on the maximal vertex 
degree. There are families of graphs where both of these magnitudes are unbounded (see Read and 
Royle [RR91], Shrock and Tsai [ST97a, 98], Thomassen [ThoOO], Bielak [BieOl], Brown, Hickman, 
Sokal and Wagner [BHSWOl], and Sokal [Sok04]). For recent discussions of some relevant research 
directions concerning zeros of chromatic polynomials and properties of their accumulation sets in 
the complex q plane, as well as approximation methods, see, e.g., Shrock and Tsai [ST97b], Shrock 
[ShrOl], Sokal [SokOla, 01b], Chang and Shrock [CSOlb], Chang, Jacobsen, Salas, and Shrock 
[CJSS04], Choe, Oxley, Sokal, and Wagner [COSW04], Dong [Don04], Dong and Koh [DK04], and 
most recently Royle [Roya, b]. 

This study of the complex roots of chromatic polynomials extends previous work that tradi- 
tionally focused on real roots, in particular, positive integer roots q which correspond to a graph 
not being properly colorable with q colors. 



6 Computational complexity connections 

The g-state Potts model partition function in Definition 2.2 involves a sum over all possible states 
of G. If G has n vertices, then there are an exponential number, g", of states. This immediately 
leads to the question of its computability. While realizing the Potts model partition function as 
an evaluation of the Tutte polynomial does not make it any easier to compute, it does enable the 
theory of one to inform the theory of the other and vice versa. In fact, it was the computational 
complexity of the Tutte polynomial in general that showed rigorously that the Potts model partition 
function is likewise intractable. The interplay of computational complexity results between the 
Tutte polynomial and the Potts model particularly illustrates the synergy between the two fields. 

We first recall the basic notions of computational complexity. A decision problem is one for 
which there is a yes or no answer, such as, can graph G be colored using k colors? P is the set 
of decision problems for which we can determine the answer in polynomial time in the size of the 
input, and NP is the set of decision problems for which we can determine if a given answer is 
correct in polynomial time in the size of the input. Whether or not P = NP remains a famous 
open question, but there is a large class of A^P-Hard problems for which finding a polynomial time 
algorithm for any one of them would automatically lead to polynomial time algorithms for all of 
them. In practice, these A^P-Hard problems are viewed to be intractable. The set of A^P-Complete 
problems are those decision problems in NP that are known to be A^P-Hard. Analogously, the set 
of #P-Complete problems is a complexity class consisting of counting problems (such as how many 
ways can a graph G be colored using k colors) that are similarly considered intractable. 

Computational complexity results for the Potts model and the Tutte polynomial have built 
in alternation upon one another as the theory has evolved. The 1990 paper of mathematicians 
Jaeger, Vertigan, and Welsh [JVW90] played a major role in this evolution. In it, the authors note 
that the Ising model {q = 2 Potts model) partition function can be reformulated as a tractable 
problem for planar graphs (referencing physicists Fisher [Fis66] and Kastelyn [Kas67]), but that 
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it was shown by Jerrum [Jer87] to be ^P-Complete in general. (Also see Vertigan [Ver05].) The 
approach in [JVW90] focuses on the problem of evaluating the Tutte polynomial along hyperbolas 
of the form {x — l){y — 1) = q where q is any real number. Note that if we let x = {q + v)/v 
and y = V + 1, and compare to Theorem 3.2, then these hyperbolas correspond precisely to the 
Potts model partition function for fixed positive integer values of q. The conclusion of [JVW90] 
is that computing the Tutte polynomial is ^P-Complete for general graphs, except when q = 1 
(which is trivial when viewed either in terms of the Potts model partition function or the Tutte 
polynomial), or when q = 2 as discussed above, or for 9 special points, namely (1,1), (0,0), (-1,0), 
{—i,i), (exp(27ri/3),exp(47ri/3)) and their reflections about the line y = x, for which the Tutte 
polynomial has easy enumerative interpretations. 

Since the Tutte polynomial, and hence the Potts model partition function, is thus typically 
computationally intractable for arbitrary graphs and argument values, a natural question arises as 
to how well either might be approximated. The answer is that in general approximating is provably 
difficult as well, but here again there is remarkable synergy between physics and mathematics, with 
results alternating between the fields. We refer the reader to an excellent overview given by Welsh 
and Merino in Section VIII of [WMOO], and to Alon, Frieze, and Welsh [AFW94, 95] for a more 
optimistic prognosis in the case of dense graphs. 

There has also been an increasing body of work since the seminal results of Robertson and 
Seymour [RS83, 84, 86] addressing computational complexity questions for graphs with bounded 
tree-width (see Bodlaender's accessible introduction in [Bod93]). A powerful aspect of this work is 
that many A^P-Hard problems become tractable for graphs of bounded tree-width. Recent research 
includes a number of results both for the classical Tutte polynomial and also for the colored Tutte 
polynomial which encompasses the Potts model with variable interaction energies. For example. 
Noble [Nob98] has shown that the Tutte polynomial may be computed in polynomial time (in fact 
requires only a linear number of multiplications and additions) for rational points on graphs with 
bounded tree width, and Makowsky [Mak05] and Traldi [Tra06] have extended this result to the 
colored Tutte polynomial. Gimenez, Hlineny and Noy [GHN06] and Makowsky, Rotics, Averbouch 
and Godlin [MRAG06] provide similar results for bounded clique-width (a notion with significant 
computational complexity consequences analogous to those for bounded tree-width; see Oum and 
Seymour [OS06]). 

While these computational complexity results for bounded tree- and clique-width are helpful in 
many instances, computing limits of the Potts model partition function as the number of vertices 
increases in an unbounded family of graphs remains an open question. Even on families of lattices 
this may be problematic. For example, Vertigan and Welsh [VW92] have shown that the Tutte 
polynomial is intractable away from q = 2 even on planar bipartite graphs (except for certain trivial 
cases such as g = or 1), and Farr [Far06] shows that computing the number of colorings of induced 
subgraphs of even the square lattice is #i-*-Complete. 

Various approaches are used to circumvent this obstacle. They include the Taylor series expan- 
sions previously discussed that provide powerful means of obtaining approximate information about 
the Potts model, as well as the computer simulations of the next section. Additionally, calculation 
of the chromatic and Tutte polynomials may sometimes be achieved for a carefully chosen family of 
graphs where the iterative operation of the deletion-contraction property leads to a solvable closed 
set of linear equations. Roughly speaking, the {m + 1) member of such a family is constructed 
by gluing a particular subgraph to the m member (see Biggs, Damerill and Sands [BDS72]). An 
example is a strip of a regular lattice of fixed width and variable length m. The resulting Potts 
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model partition function has the form of a finite sum of m powers of a set of algebraic functions 
multiplied by certain coefficients. These algebraic functions are the roots of a set of equations 
resulting from the iterative operation of the deletion-contraction theorem or equivalently, are the 
eigenvalues of a certain type of transfer matrix. Some calculations of chromatic polynomials of re- 
cursive families of graphs include Biggs, Damerill and Sands [BDS72], Beraha and Kahane [BK79], 
Beraha, Kahane, and Weiss [BKW80], Read [Rea88], Salas [Sal90, 91] , Read and Royle [RR91], 
Shrock and Tsai [ST97a], Rocek, Shrock and Tsai [RST98], Shrock and Tsai [ST99], Shrock [Shr99, 
01], Biggs and Shrock [BS99], Sokal [SokOO], Salas and Sokal [SSOl], Chang and Shrock [CSOla, b], 
Biggs [BigOl, 02a], Jacobsen, Salas and Sokal [JSS03], Jacobsen and Salas [JSOl, 06]. 



7 Monte Carlo simulations of the Potts model 

The computational intractability of the Potts model partition function has led to the development of 
Monte Carlo simulations for the model; see the texts of Newman and Barkema [NB99] and Landau 
and Binder [LBOO] for additional background on the methods described below. We illustrate the 
basic principles of this kind of simulation in the simplest, q = 2, case of the Ising model on a square 
lattice, and then briefly mention some modifications leading to more sophisticated simulations. 

Since complexes are often very large, with many different spin choices for their elements, the 
probability of a single state appearing out of the exponential number of states is nearly zero, but 
the macroscopic properties for many different states may be similar. Therefore, the goal is to 
determine the average characteristics the system is likely to exhibit in the long run; i.e., we want 
to approximate the expected value of a macroscopic property when the system is in equilibrium. 
In the case of the Ising model, we might want to determine the expected value of the magnetism 
at a given temperature. The simulation must compute this expected value by averaging over a 
sufficiently large sample of states that correspond to an independent random sample of states from 
the Boltzmann distribution. These states are generated through a Markovian random walk on the 
lattice. 

It is preferable, but not necessary, to begin with an initial state that is characteristic of the 
temperature at which the properties of the system are being measured. For example, if one were to 
start with an ordered spin configuration at a high temperature, then considerable computer time 
would be expended to warm up the simulation, while if one starts with a random spin configuration, 
much less time is spent reaching equilibration. When the system is at equilibrium, the value of 
the macroscopic property of interest should stay within a fairly small range. The simulation is 
generally run from a number of different initial configurations to ensure that the system has actually 
found the equilibrium value, rather than a locally stable value. Since a state is clearly dependent 
on a few of the previous states in the random walk, an autocorrelation function is computed to 
determine the distance between samples taken in the random walk to ensure that the sample 
points are independent. The necessary simulation length and corresponding statistical error can 
then be estimated in the typical manner for applications of the Central Limit Theorem. As with 
all experiments, systematic error may occur and can be difficult to detect. 

Prom the initial state, each vertex is visited in turn, and the program computes the probability 
ratio comparing the likelihood of the vertex changing its spin versus retaining its current spin. This 
simulation captures the effect of temperature on the model, encoding the tendency of the system to 
move toward a lower energy state at low temperatures and remain agitated at high temperatures, 
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as follows. Recall that the probability of a state occurring is Pr (w) = (exp {—[3hi {w)) /Pi (G), so 
the ratio of the probability of a new state Sm to the probability of the current state Sc is 

exp(-/3/ii(5Ar)) 



Vt{Sn) _ J]exp(-/3/ii(u;)) _ ex.v {- (ihi{S n)) _ ( hi{Sc) - hi{SN , 



Vt{Sc) exp(-/?/i,(gc)) exp(-/?/i,(5c)) ^V ^^T 

X;exp(-/3i(a;)) 

Note that this avoids computing the generally A'^P-Hard partition function, Pj. Also note that 
since Sn differs from Sc only in a change of spin at one vertex u, the computation of hi (Sn) is 
exactly the same as that for hi (Sc) at every edge except for those incident with v. 

In the commonly used Metropolis Monte Carlo algorithm, if the new state has lower energy 
than the current state, hi (Sn) < hi (Sc), the algorithm changes the system from state Sc to state 
Sn- However, if hi (Sc) < hi (Sn), the program compares 

Ft{Sn) (hi{Sc)-hi{SN) 

— exp ' 



Pr(5c) V i^T 

to a random number r with < r < 1 and changes state if r < Pr(S'/vr)/Pr(S'c). At high 
temperatures, this ratio will be nearly 1 regardless of the Hamiltonians, so spins will continue 
changing with negligible preference for lower energy states. On the other hand, if the temperature 
is quite low, the system strongly favors low energy states. 

Although the behavior of the model is clear at very high or very low temperatures, it is less 
apparent what happens at midrange temperatures. One of the fundamental questions for the 
Potts model on a regular d-dimensional lattice is determining the critical temperature Tc for 
a phase transition. In Monte Carlo simulations, temperatures near the critical value can cause 
computational challenges, due to the increased statistical error in that range as well as the increase 
in the autocorrelation time. Nevertheless, simulations are an important tool in the study of the 
Potts model, as few exact analytic results are known. 

There are various refinements of this basic model leading to more sophisticated simulations. 
For example, there is no need for the underlying graph to be a square lattice. It can be a different 
regular lattice, such as triangular or honeycomb in two dimensions; a cubic, face-centered cubic, 
or body-centered cubic in three dimensions. More generally, it can be any graph appropriate to 
the application, even a complete graph if every site interacts with every other, although dense 
and/or irregular graphs can present programming challenges. The simulation can be extended to 
larger numbers of spins by computing the relative probabilities comparing the current spin at a 
vertex with each of the other possible spins. In the heat-bath algorithm, the probability ratios 
are normalized so they sum to 1, and then each is assigned a proportional segment of the unit 
interval. A random number is generated in the unit interval, and the spin is changed according 
to the segment that contains the random number. There are also useful techniques for improving 
the speed of these simulations, including clustering methods. Clusters of locally aligned states slow 
down the simulation, since the likelihood of a flip occurring is very low. Therefore, the simulation 
will spend extended periods in the same state. The Wolff algorithm improves the running time by 
flipping these clusters of like spins together instead of considering them one by one. See Wang, 
Kozen, and Swendsen [WKS02], Deng, Garoni, and Sokal [DGS07], and Deng, Garoni, Machta, 
Ossola, Polin, and Sokal [D-l-07]. 
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Other significant modifications and variations of the model include inclusion of an external mag- 
netic field, next-nearest-neighbor interactions, edge-dependent interaction energies, and additional 
terms in the Hamiltonian. In realistic physics studies, one also must often consider the effects of 
disorder in actual materials, such as vacancies and impurities in a crystal lattice. We give further 
discussion of some of these variations in the next section. 



8 Why is this model attracting so much attention? 

Besides its intrinsic mathematical interest, the Potts model, in many variations, is increasingly 
applicable to a wide variety of complex systems where local interactions can predict global behavior. 
This is particularly true as computing power has enabled increasingly powerful and predictive 
simulations and as researchers have found sophisticated modifications of the model to more closely 
mimic the behaviors of various systems. The popularity of the Potts model is roughly indicated by 
a recent Google Scholar search for "Potts model" producing over 63,000 hits. We give a sample of 
applications here, just to demonstrate the scope of this theory. In these examples, the Hamiltonian 
is extended to encode forces in addition to simple nearest-neighbor interactions, but the probability 
distribution, and hence partition function, is still generally defined analogously to that of Definition 
2.2. 

The original magnetism application addressed by Ising considers the overall ferromagnetic ( "nor- 
mal" magnetism) behavior of a lattice where the two possible spins at each position are positive 
and negative. The energy of the system is minimized if all points on the lattice have the same 
spin, while, to maximize entropy, all states should be equally likely (which would strongly favor 
nonmagnetic states). The Boltzmann distribution quantifies the relative importance of energy and 
entropy in determining the likelihood of a given state in terms of temperature (assuming the system 
is at the same temperature as its surrounding environment). The standard Hamiltonian is given in 
Definition 2.1, but is sometimes extended to include an external magnetic field: 

h{uj) = -J ^ 5{ai,aj) - H'^ai. 

ij€E(G) i 

Simulations and series expansions are used in higher dimensions to determine the phase transi- 
tion temperature below which the system exhibits a nonzero spontaneous magnetization and above 
which this magnetization vanishes. See Stanley [Sta71], McCoy and Wu [MW73], Chandler [Cha87], 
Plischke and Bergesen [PB06] for the Ising model, and Stanley [Sta71], Carter [CarOl], Plischke 
and Bergesen [PB06] for further information on the Boltzmann distribution. 

Sanyal and Glazier [SG06] employ the Potts model to simulate foam flow and investigate in- 
stabilities, i.e., velocities at which larger bubbles start to flow faster than smaller bubbles. Two 
adjacent lattice points have the same spin if and only if they are part of the same bubble, hence 
the number of spins Q is extremely large. The Hamiltonian is: 

Q 

ijeE{G) n=l 

where J is the coupling strength at the boundary between two bubbles, an is the actual area of 
the bubble. An is the area the bubble would have if it were not subjected to external forces, and A 
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is the strength of the area constraint on the bubble (based on the compressibiHty of the gas). At 
each step in the simulation, the spin at a randomly selected lattice site is considered; if that site 
is along the boundary with another bubble, a switch in spin of that lattice site to the neighboring 
bubble is considered (and accepted with some probability). 

Turner and Sherratt [TS02] use an extension of the Potts model to study cellular malignancy 
growth. They are particularly interested in the impact of the relative strength of a few factors 
known to affect cell behavior. The Hamiltonian for their model is: 



Here ij indexes the lattice point in the i*'^ row and j*'^ column of a two dimensional grid. The 
spin value aij records which of the many cells contains the ij lattice point. The first term of the 
Hamiltonian encodes at the interaction energy between a cell and its eight nearest neighbors on 
the lattice. Adjacent lattice sites with the same spin value aij represent a single cell, so there will 
be no interaction in that case. Otherwise the interaction energy depends on the cell type r^aij), 
which may be normal or malignant. The second term models the energy required for a cell a to 
maintain a volume different from its natural volume in the absence of external forces, similar to the 
foam example above. Instead of temperature, the /3 in this application corresponds to a diffusion 
coefficient affecting the random motility of the cells. Turner and Sherratt [TS02] further extend 
the Potts model by allowing cells to replicate, hence changing the lattice, during the simulation. 
In doing so, they take into account the interaction of a cell with the extracellular (protein) matrix. 

Nobel laureate Thomas Schelling published a seminal paper titled Dynamic models of segregation 
in 1971 [Sch71] that considers the possibility of micro-motive explanations for racial segregation 
(in addition to organized and economic explanations). The premise is that individual decisions to 
avoid minority status (or to require being in a minority of some minimum size) could lead to the 
macro-effect of segregation. Schelling places vacancies, stars, and zeros randomly on a checkerboard 
and then iteratively considers the happiness of the stars and zeros with their local neighborhoods, 
moving an unhappy star or zero to the nearest vacant spot that meets their happiness criteria. 
Meyers-Ortmanns [Mey03] models a similar premise to Schellings (that micro-motive explanations 
can lead to immigrant ghettos) with a more Potts-like model where the Hamiltonian measures the 
happiness of individuals with their neighbors, the temperature is viewed as a social temperature 
where warmer temperatures reflect facilitation of integration and assimilation, and at each step in 
the simulation two neighbors are able to exchange places with a probability based on the likelihood 
of the new state with respect to the current state. Schulze [Sch05] extends Meyers-Ortmanns to 
address up to seven different ethnic groups. 



9 Some active areas at the interface of combinatorics and statis- 
tical mechanics 

While there is intense interest in the Potts model, not just from the physics community, but 
throughout the sciences, its properties are of intrinsic combinatorial interest as well, and clearly 
much work remains to be done in exploring, developing, and extending them. For readers who 
would like to pursue broad perspectives and further background on the topics mentioned here, 
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overviews of the Tutte polynomial may be found in Brylawski [BrySO], Tutte [Tut84], Brylawski 
and Oxley [B092], Biggs [Big93], Welsh [Wel93, 99], Bollobas [Bol98], and Farr [Far07], with the 
latter three also discussing the Potts model from a mathematical point of view. Some works 
containing reviews of the Potts model from physics and mathematical physics viewpoints include 
those of Baxter [Bax82], Wu [Wu82, 88], Cipra [Cip87], Martin [Mar91], Welsh and Merino [WMOO], 
Shrock [ShrOl], Chang, Jacobsen, Salas, and Shrock [CJSS04], Sokal [SokOO, 05], and Farr [Far07a]. 

In addition to the areas already discussed, researchers are also currently very interested in a 
number of related areas. New computational techniques for relevant combinatorial polynomials, 
particularly those that can be applied to regular structures such as lattices, are always sought. Of 
independent, yet mutual, interest for both graph theory and statistical mechanics is the determi- 
nation and physical interpretation of the zeros of Tutte polynomial and its cognates. The effects of 
boundary conditions on lattices, including periodicity leading to toroidal and other topologies, are 
an important area of study, as are the connections to knot theory. Improved Markov-chain Monte 
Carlo methods are in constant demand, as are further computational complexity results. Recent 
work also includes the use of methods from statistical mechanics in combinatorial enumeration. 
Space prohibits providing an exhaustive list, but a few examples not previously mentioned in these 
areas include Chang and Shrock [CSOlc]; Sokal [SokOla]; Wu and Wang [WWOl]; Woodall [Woo02]; 
Traldi [Tra02, 06]; Bonin and de Meir [BdMN03]; Jackson [JacOS]; Read [Rea03]; Biggs [Big04], 
Royle and Sokal [RS04]; Morris and Peres [MP05]; Makowsky [Mak05]; Wu [Wu05, 06]; Makowsky, 
Rotics, Averbouch, and Godhn [MRAG06]; Oum and Seymour [OS06]; Farr [Far07b]; Martinelh, 
Sinclair, and Weitz [MSW07]; Jerrum [Jer07], and Jacobsen and Salas [JS07]. 

We close with an observation that, although far from original, hopefully may encourage com- 
binatorialists to reach across the dark spaces between disciplines and engage in research related to 
statistical mechanics. Theoretical physicists have produced a wealth of information about phase 
transitions and critical phenomena leading to well-supported assertions, many of which still need 
rigorous mathematical treatment or lead to questions of intrinsic mathematical interest. The result 
is a ready supply of appealing mathematical problems. This is especially true for combinatorialists, 
into whose domain many of these problems from statistical mechanics naturally fall. 
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